*This file performs the bootstrap analysis referenced in Footnote 18. 
*The results are plotted in Appendix Figure A1.

********************************************************************************
*DEFINE DIRECTORIES
local home CHILD
local main CHILD/JPE
local logs CHILD/JPE/logs
local data CHILD/JPE/data
local results CHILD/JPE/results
local network NETWORK
********************************************************************************

cd "`main'"

******Section A: Replicate sample pull
*generate log file
cap log close
local time = subinstr(c(current_date)," ","",.)
log using bootstrap_outcomes.txt, text replace

*bring in data
cd "`main'"
use regression_data_final_jpe.dta, clear

*merge in single treatment variable
preserve
use all_summary_statistics.dta, clear
keep mmi_id treated-badfda nd adhd
tempfile temp
save `temp',replace
restore

drop treat
drop drug
drop bad

merge 1:1 mmi_id using `temp'
keep if _merge==3
drop _merge

*get rid of extra variables
drop tca benzo badfda

*define drug treatment outcome
generate drugs = drug_only+both>0
drop drug_only therapy_only both

*merge in new bad prescribing variables
merge 1:1 mmi_id using bad_prescribing.dta
foreach var of varlist fda_ok-red_flag {
	replace `var' = 0 if `var'==.
}
drop _merge

*reorder variables
order treated drugs red_flag grey_area fda_ok, after(mmi_id)

*merge in REVISED DIAGNOSIS INFO
cd "`home'"
merge 1:1 mmi_id using diagnosis_revised.dta
replace ada = 0 if _merge==1
replace exclude = 0 if _merge==1
generate mood_adj_anx = ada
drop _merge

*POPULATION SELECTION
keep if mood_adj_anx==1 & exclude==0
drop exclude ada mood_adj_anx

*require at least 10 kids in a zip3 zone
by zip3, s: generate count2 = _N
drop if count2<10

*parse variables
keep mmi_id-fda_ok zip3

*save file at this point
tempfile data
save `data', replace

******Section B: Identify into groups (see iPad code)

*define the treatment bin
generate group = ""
replace group = "A" if treated==0
replace group = "B" if treated==1 & drugs==0
replace group = "C" if red_flag==1
replace group = "D" if grey_area==1
replace group = "E" if fda_ok==1

*parse variables
keep mmi_id zip3 group

*save this file
tempfile temp
save `temp', replace

******Section C: Run the bootstrap simulation
*note, with strata, you get the right N within each strata BUT you are sampling WITHIN the strata
*this version of Section C uses the packaged Chi2 test in Stata---much, much faster!

*create the backbone
use `temp', clear
keep zip3
generate counter = _n
tempfile backbone
save `backbone', replace

*run the bootstrap
use `temp', clear
drop mmi_id
quietly tab zip3 group, chi2
generate chi2 = `r(chi2)'
keep if _n==1
drop zip3 group 
generate iteration = 1
tempfile temp2
save `temp2', replace

forvalues i = 2(1)10000 {
use `temp', clear
keep group
bsample
generate id = runiform()
sort id
drop id
generate counter = _n
quietly merge 1:1 counter using `backbone'
drop _merge
quietly tab zip3 group, chi2
generate chi2 = `r(chi2)'
quietly keep if _n==1
drop zip3 group counter
generate iteration = `i'
append using `temp2'
quietly save `temp2', replace
di "`i'"
}
cd "`main'"
save bootstrap_test_fast.dta, replace

*draw graph that is in the paper
histogram chi2, lcolor(none) color("158 202 225") width(5) xline(999 2947,lpattern(dash) lcolor("49 130 189")) xlabel(0 500 999 `" "Critical value" "(999)" "' 1500 2000 2500 2947 `" "Test statistic" "(2947)" "', labsize(vsmall))  xtitle("Chi-square statistic bootstrap", size(small)) title("Empirical Distribution of Bootstrap Chi-square Values", size(medium)) note("N=10,000 bootstrap replications." "Degrees of freedom = 928.", size(tiny)) graphregion(fcolor(white) lcolor(none)) ytitle(,size(small)) ylabel(,labsize(vsmall))
cd "`results'\Final Tables"
graph export bootstrap.pdf, replace as(pdf)